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Abstract 

Reaction-diffusion problems are often described at a macroscopic scale 
by partial derivative equations of the type of the Fisher or Kolmogorov- 
Petrovsky-Piscounov equation. These equations have a continuous family 
of front solutions, each of them corresponding to a different velocity of the 
front. By simulating systems of size up to N = 10 16 particles at the micro- 
scopic scale, where particles react and diffuse according to some stochastic 
rules, we show that a single velocity is selected for the front. This veloc- 
ity converges logarithmically to the solution of the F-KPP equation with 
minimal velocity when the number N of particles increases. A simple 
calculation of the effect introduced by the cutoff due to the microscopic 
scale allows one to understand the origin of the logarithmic correction. 

PACS numbers: 02.50.Ey, 03.40.Kf, 47.20. Ky 



1 The Fisher equation 

The Fisher equation [Q, also called KPP equation (for Kolmogorov-Petrovsky- 
Piscounov^J) is widely used to describe front propagation in many problems of 
physics, chemistry and biology: 

dc d 2 c 2 . . 

Fisher first introduced this equation to represent "The Wave of Advance of 
Advantageous Genes" fl in a population. The concentration c(x, t) was the 
fraction of individuals in a population at position x and time t that exhibit 
some benefic genes, and (Q) was used to describe how this favorable gene would 
spread in the population. Equation ([!]) can also model the dynamics of sick 
individuals in a population during a viral contagious infection, the proportion of 
burnt-out gases in a combustion J3|, the concentration of some species produced 
in a chemical reaction, etc. It also appears in the mean field theory of directed 



1 



polymers in a random medium B and in the calculation of Lyapunov exponents 
of large sparse random matrices || ||. 

In (|l|), c(x,t) represents a concentration, implying that, for all x and t: 

0<c(x,t)<l. (2) 

If we look at solutions constant in space (dc/dx — 0), equation ([!]) becomes 
simply: 

There are two stationnary solutions: c — (unstable fixed point) and c = 1 
(stable fixed point). 

The role of the diffusion term d 2 c/dx 2 in ([[]) is to spread any positive pertur- 
bation. Therefore, if initially c(x, 0) > in some region of space and c(x, 0) = 
elsewhere, the perturbation will grow, reach asymptotically the stable fixed 
point c = 1 and spread throughout the whole space. Ultimately, the stable 
value c = 1 will be reached everywhere. 

To study how the stable region (where c = 1) invades the unstable one 
(where c = 0), one can consider an initial condition where c(x,0) decreases 
monotonically from c(— oo, 0) = 1 to c(+oo,0) = 0. If at time t = the initial 
condition c(x,0) decreases fast enough as x — > +oo (in particular if c(x,0) is 
a step function), the front moves in the long time limit with a well defined 
speed v min [@: 

^min = 2. (4) 

To understand from (||) why the velocity v m i n — 2 is selected looks a priori 
a very hard task. Equation (Q) is a non-linear partial derivative equation and 
there is no way of writing the full expression of c(x, t) for an arbitrary initial 
condition. However, the velocity w m i n = 2 can be understood easily without 
solving the full non-linear problem: let us assume that the front moves at some 
constant speed v. The concentration profile c{x, t) takes then the form: 

c(x, t) = F v (x — vt), (5) 

where F v (z) is the solution of: 

F" + vFl +F V - Fl = 0, (6) 

with F v {— oo) = 1 and F v (+oo) = 0. In the region where F v {z) 1, i.e. far 
ahead of the front, one can neglect in (^|) the non-linear term. The general 
solution of the linearized equation is a sum of two exponentials so that for z — > 
+oo one of these two exponentials dominates: 

F v (z) ~ Ae-r, (7) 

and, from the linearized version of (^|), one finds that v and 7 are related by: 

«( 7 )= 7 +I. (8) 

7 
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We see that the asymptotic decay 7 in ^ determines completely the velocity 
^(7) of the front. 

What (||) tells us is that any velocity v > u m i n (u m i n = 2) is possible for the 
front. (It should be noted that v < v m - ln would also be possible by allowing 7 
to be complex. However, a front moving at such a speed would take negative 
values in the tail, and this would violate condition (|2|).) 

The minimal speed u m i n = 2 reached for 7, n i n = 1 has a special status: it 
has been shown ||, |[ |Io|, [n], [l2| that if in its initial condition the front decays 
faster than e~ 7miIlX (in particular, if c(x,0) is a step function), then the front 
moves asymptotically with this minimal speed w m i n . Moreover, in the long time 
limit, the position X(t) of the front is given by 

X(t) = 2t- -hit + 0(1). (9) 

(In other words, the velocity of the front converges to v m i n — 2 with a leading 
correction fll2^ given by — 3/(2i). The presence of a logarithmic correction in (||) 
makes often a precise determination of the asymptotic velocity difficult.) 

Most properties of ([!]) (selection of the minimal velocity for steep enough 
initial conditions, logarithmic corrections to the position as in (^)) can also be 
recovered in a whole class [[lO] of front equations where a stable region invades 
an unstable one. An example very different from (^) that we will consider below 
is: 

c(x, t + 1) = 1 - [1 - J da p{a)c(x - a, t)] 2 , (10) 

where p(a) can be any density function (p(a) > and J da p(a) = 1). As 
for (|J), the uniform solutions c = and c = 1 are respectively unstable and 
stable and the integral over a spreads any positive perturbation as does the 
diffusion term in (Q). 

As for (|l|), the linearized version of ([h]) where terms quadratic in c are 
neglected determines the velocity. For an exponential decay (]?]) of the front, 
one finds a dispersion relation 1^(7) generalizing (0): 



vM = — In 

7 



2 / da p(a)e< 



(11) 



and for a steep enough initial condition the minimum velocity 



Vmin = minw(7) = w(7 min ) (12) 

7 

is reached in the long time limit. The position X(t) is then given Jll| for large t 
by: 

X(t)=v min t--^—\nt + 0(l). (13) 

Tin in 

(Note that in general 7 m j n in (|l2|) is finite except for very particular choices of 
p{a).) 
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2 The microscopic stochastic model 



Front equations of type (|l|) or (|T^) originate often as the large-scale limit of 
microscopic stochastic models || ^, |l3|, [FI], [TJ|. Here we study a particular 
microscopic model which, as we will see, is described in the large scale limit by 
the front equation (fTo|). We will compare the velocity measured for this micro- 
scopic stochastic problem with the velocity ( p| , ^2|) expected for the traveling 
wave equation (fio|). 

Our microscopic model [jlT|| is defined as follows: Imagine a population where 
each generation has exactly N individuals. Each individual i (1 < i < N) at 
generation t (t is an integer) is characterized by its fitness Xi (t) , a real number 
representing its adaptation to the environment. The state of the system at any 
time t is completely determined by the N numbers Xj(i). 

At time t = 0, we set 2^(0) = for all i (but this choice of initial condition 
is actually unimportant in the long time limit) . By definition of the model, the 
Xi(t) evolve from generation t to generation t + 1 with the following rule: 

Xi(t + 1) = max [x mi (t) + a i} Xf i (t) + a[] , (14) 

where m, and fi are the two parents of the new individual i, chosen at random 
in the previous generation t (in other words, rrii and fi are random integers 
uniformly distributed between 1 and N), and where and are random 
numbers independently chosen according to some probability distribution p(a) 
representing random mutations. So at each generation, the rrii, fi, ai and a[ 
are independent and new values are chosen at every time step. 

Under the dynamics (|l4|) , the cloud of N points Xi (t) moves along the line 
and we want to determine its asymptotic velocity v^, that is: 

v N = hm — — , (15) 

t^+oo t 

where 

1 N 

Let us now see how one can relate this microscopic model to the traveling 
wave equation (f[o|). We define c(x, t) as the fraction of population which has a 
fitness larger than x: 



1 - 

<x,t) = -j2 @ ( x i( t )- x )- ( 17 ) 



(By convention, we choose here Q(x) = 1 if x > and 0(x) = if x < 0.) 
Obviously, c(x, t) is a monotonic decreasing function of x going from c(— oo, t) — 
1 to c(+oo,<) = 0. At time t = 0, we have c(x,0) — 1 for x < and c(x,0) = 
for x > 0, so the initial condition is a step function. It should be noted 
that c(x,t) can only take values which are integral multiples of l/N. 
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Clearly from (|T7j), the position X(t) of the cloud of points can be rewritten 
with the function c(x, t) as: 

/+oo 
dx [c(x,t) - c(x,0)] . (18) 
-OO 

Given the positions Xi(t) of all the particles (or, equivalently given the func- 
tion c(x, t)), the Xi(t + 1) obtained from ( |l4| ) are independent random variables. 
Therefore, if we fix c(x, t), the average (c(x, t + 1)^ over the dynamics between 
time t and t + 1 gives: 

<c(z, t + l>> = l-[l-/d Q p(„)c(*- <,,*)]*, (19) 

which, except for the ( ), is exactly ([!(]). However, if we try to average over the 
whole history (i.e. over all the timesteps), we need to average terms quadratic 
in c on the right hand side of (|l9|). This means (c(x,t+l)) is not only related to 
(c(x, t)), but also to correlations like (c(x, t)c(x',i)), and this makes the problem 
very difficult for finite N. On the other hand, if we neglect these correlations 
(and one can argue that these correlations are small for large enough N) and 
replace (c(x, t)c(x' , t)) with (c(x, t))(c(x', t)), then ( |l9| ) reduces exactly to (|l0|). 
So ( |l0| ) can be thought as the mean-field (or large N) version of the microscopic 
model Q. 

As the initial condition c(x, 0) given by ( |l7| ) is a step function, the mean field 
equation predicts a front moving at the minimum velocity v^in given by ( O |l2| ). 

3 Direct simulations 

We have simulated the microscopic model (|l4]) for several choices of the distri- 
bution p(a): the uniform distribution in the range < a < 1 

Puni(a) = e(a)0(l - a), (20) 
the exponential distribution 

Pcxp(a) = 6(a)e- Q , (21) 

and a discrete distribution 

Wisc (a) = 0.25 (5(1 - a) + 0.75 5(a). (22) 

We have also simulated a generalization of the problem (Martian genetics) where 
each new individual at time t + 1 has three parents, so that (Q) is replaced by 
the max over three terms, with the effect of mutations distributed according 
to ©. 

The minimal value of the speed w m i n and the corresponding decay rate 7 m i n 
of the deterministic front equations ([h], [ll], |l2|) for these four cases are given in 
table 0. 
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model 


P = Puni 


P = Pcxp 


P = Pdisc 


Martian model 


Tin in 


5.262 076... 


0.626 635.. . 


2.553 245.. . 


8.133 004... 




0.815172... 


2.678 347... 


0.810710... 


0.877338.. . 



Table 1: Values of 7 m i n and v m - nl for different models. 



For these four models, we have simulated (|l4j) for T = 10 7 timesteps after a 
transient time of T' = 10 6 timesteps to eliminate the effect of initial conditions. 
For several choices of N, we have measured the speed as: 

X(T + T)-X(T) 
V N = y • (23) 

Figure [l] represents the difference between the mean-field speed i> m in (given in 
table |l|) and the speed vn measured in the simulation for several choices of N 
(16, 32, 64, 128, 256, 512, 1024, 2048 and 4096). 




Figure 1: Log-log plot of the difference between the speed v m - m given by the 
mean-field theory and the speed vn measured in Monte-Carlo simulations for 
the four models as a function of the number N of particles. The straight line 
represents iV -1 / 3 . 

We see on figure |l| that a single speed is selected in the microscopic 
model. This speed is always lower than the speed v m - m of the deterministic 
front equation, and the difference v m \ n — seems to decay like iV -1 / 3 for all 
the variants of the model. The effective power law exponent seems however to 
decrease slowly as N increases. 

In order to confirm the TV -1 / 3 decay of figure (|l|), we tried to increase AT, 
but as in the direct simulation computer time scales like N x T, it is very hard 
to make N much larger than 10 4 or 10 5 for a number T = 10 7 of timesteps. 
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4 Highly parallel simulations 

We are now going to describe a computational trick which we developped |TT[ 
for some particular distributions p{a) such as 

p(a)=p6(a-l) + (l-p)6(a), (24) 

allowing to simulate the microscopic model for a huge number of points (N up 
to 10 16 ). We restrict p to the range < p < 0.5 (to avoid one of the rare cases 



where (11) has no minimum for a finite 7 m i n )- 

For the distribution (pi]), the Xi(t) are always integers if they are so at t = 
and the concentration c(x, t) as defined by (|l7| ) is constant between any pair of 
consecutive integers. We call respectively x m i n and x max the positions of the 
leftmost and rightmost particles at time t and w = x max — x m j n + 1 the width 
of the front. We observed in our simulations that w is typically of order In TV, 
so that even for N as huge as 10 16 , the number of possible values of the Xi{t) 
at a given time is very limited, and the whole information in c(x, t) is carried 
by the number of particles at each integer x between a; m ; n and x max . 

Knowing the function c(x, t) at time i, we generate c(x, t+l). As at time t, all 
the Xi(t) satisfy £ m in < Xi(t) < x max , from ( fl4| ) and ( ]24| ) the positions Xi(t + 1) 
will lie between x m in and x max + 1 • The probability pk that a given particle i 
will be located at position x m i n + k at time t + 1 is 

Pk = (c(x min + k - l,t + 1)> - (c(x min +k,t+ 1)), (25) 

with (c(x,t+ 1)) given by (|l9|). Obviously pk ^ only for < k < w. The 
probability to have, for every k, particles at location x m ; n + k at time t + l 
is given by 

TV' 

P(n ,n u ...,n w )= ; P%° pT ■ ■ ■ P? w 

"iln, I 7J....I 



no! nil 

xS(N — no — ni — n w ). (26) 



Using a random number generator for a binomial distribution|16| , expression (|2 
allows us to generate the random numbers nt directly fll)] . 

We have measured from (|2^) the velocity vn of the front for several choices 
of p (0.05, 0.25 and 0.45) and for N ranging from 100 to 10 16 . Figure || shows 
the results together with functions A(p) j In 2 N given in the next section. 

Clearly the apparent power law of figure ^ does not persist as N increases 
and the simulations indicate that — vn ~ ln~ 2 TV for large N. We are going 
to see in the next section that this logarithmic correction has a simple origin. 



5 Effect of the cutoff 



The two main differences between the traveling wave equation (10) and the 
microscopic model (14) is that the microscopic model is stochastic and that 
c(x, t) varies by steps multiple of l/N (see (p"7j)). 
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N 



Figure 2: Differences between the speed obtained by the mean- held theory and 
the speed measured in highly parallel simulations for different values of p as a 
function of the number N of particles. The lines represent for each value of p 
the prediction (p7|). 



The effect of the noise is hard to treat analytically and we have not succeeded 
yet to make a satisfactory theory of it. The effect of discretization can be 
however understood rather simply: let us modify (|l0| ) by imposing, after each 
timestep, c(x, t + 1) — if the value given by PQ ) is smaller than 1 /N, (in other 
words, we put a cutoff in the deterministic model to mimic the fact that c(x, t) 
changes by steps of 1/N). The velocity v' N of this deterministic model with a 
cutoff can be easily measured as under these dynamics the front reaches rapidly 
a periodic regime jl^]. As the cutoff goes to zero (i.e. as N goes to infinity) , 
the speed v' N converges to the mean-field speed w m i n given by ([□]) and (|l2|). 
On figure |we compare this speed for N = 64 and N = 512 with w m i n for p(a) 
given by ( p4| ) with < p < 0.5. (One can note that the speed gets locked to 
rational values as p varies.) 

The speed v' N of this new deterministic model can be calculated analyti- 
callyfl) for large N: 



T 2 7^un u " (Train) 

' J N " ' 



(27) 



21n^7V 

where u('y) is given by (|TT|) . Comparing the prediction ( j27] ) with the results of 
the simulations in hgure^ gives a good, though not perfect, agreement. This 
indicates that the slow convergence of the velocity of the stochastic model is 
controled by the effect of the cutoff. 



6 Conclusion 

We have seen that, for a very particular microscopic stochastic model (^4|), it was 
possible to simulate systems of 10 16 particles. Our model is so particular that 
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Figure 3: Comparaison of the velocity v m i n valid for infinite N (top curve) with 
the velocity of the model with a cutoff 1/N for TV = 5f2 (middle curve) and 
N = 64 (bottom curve) as a function of p. 

there is no hope that our trick could be extended to large classes of statistical 
physics problems. In our case, however, going to very large N enabled us to 
clearly discriminate between a power law and a logarithmic correction. 

The fact that the microscopic scale selects a single velocity]!!], [Tt], |l8| with a 
logarithmic correction due to a cutoff seems to appear in several related problems 
(reaction-diffusion || kinetic theory Q). Even model ( pl| ) can be introduced 
in many different contexts, like directed polymers in a random medium (where 
Xi(t) would be the free energy of a polymer of length t ending at position i), 
or growth problems (where Xi(t) would be the height variables). Still, we do 
not know yet whether the prediction (^) based simply on the effect of the 
cutoff gives the exact large- N behavior of the model (|14|), or whether a more 
sophisticated theory is needed to explain the results of section ||. 

We thank M. Droz for communicating us several relevant references. 
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